# Clean workspace
rm(list = ls())  # Remove all objects
graphics.off()    # Close all graphics devices
cat("\014")      # Clear console

# Install required packages if not already installed
if (!require(tidyverse)) install.packages("tidyverse")
if (!require(readxl)) install.packages("readxl")
if (!require(writexl)) install.packages("writexl")  # Add writexl package for Excel output
if (!require(moments)) install.packages("moments")  # Add moments package for skewness and kurtosis
if (!require(nortest)) install.packages("nortest")  # Add nortest package for normality tests

# Load required libraries
library(tidyverse)  # for data manipulation
library(readxl)     # for reading Excel files
library(writexl)    # for writing Excel files
library(moments)    # for skewness and kurtosis calculations
library(nortest)    # for Anderson-Darling normality test

# Function to load data
load_behavior_data <- function(file_path, sheet_number = 3) {  # Default to load the 3rd sheet
  # Determine file extension
  file_ext <- tools::file_ext(file_path)
  
  # Load data based on file type
  data <- if (file_ext == "csv") {
    read.csv(file_path)
  } else if (file_ext %in% c("xlsx", "xls")) {
    read_excel(file_path, sheet = sheet_number)  # Specify sheet
  } else {
    stop("Unsupported file format. Please provide a CSV or Excel file.")
  }
  
  # Convert to tibble for better handling
  data <- as_tibble(data)
  
  # Ensure proper column types
  data <- data %>%
    mutate(
      Subject = as.factor(Subject),
      Semantic_Relatedness = as.factor(`Semantic Relatedness`),
      Priming_Word_Type = as.factor(`Priming Word Type`),
      Group = as.factor(Group)
    )
  
  return(data)
}

# Set input file path
input_file <- "source_data/Ch_Ja_Lex_Behavior_2Group_29Sub_for R.xlsx"

# Create result folder name from input file name
result_folder <- tools::file_path_sans_ext(basename(input_file))
result_folder <- gsub("_for R", "", result_folder)  # Remove "_for R" suffix
result_path <- file.path("results", result_folder)

# Create result folder
if (!dir.exists(result_path)) {
  dir.create(result_path, recursive = TRUE)
  print(paste("Created result folder:", result_path))
}

# Load your data (3rd sheet)
print("Loading data from 3rd sheet...")
behavior_data <- load_behavior_data(input_file, sheet_number = 3)
print("Data loading completed")

# View data structure
print("\n============= Data Structure =============")
print("Data dimensions:")
print(dim(behavior_data))
print("\nVariable names:")
print(names(behavior_data))
print("\nData structure:")
str(behavior_data)

# View factor levels
print("\n============= Factor Levels =============")
print("Number of subjects:")
print(length(unique(behavior_data$Subject)))
print("\nSemantic Relatedness levels:")
print(levels(behavior_data$Semantic_Relatedness))
print("\nPriming Word Type levels:")
print(levels(behavior_data$Priming_Word_Type))
print("\nGroups:")
print(levels(behavior_data$Group))

# Analyze RT normality
print("\n============= RT Normality Analysis =============")

# Check if RT column exists
if ("RT" %in% names(behavior_data)) {
  print("RT data found. Analyzing normality...")
  
  # Basic descriptive statistics
  rt_summary <- summary(behavior_data$RT)
  print("RT Summary Statistics:")
  print(rt_summary)
  
  # Skewness and Kurtosis
  rt_skewness <- skewness(behavior_data$RT, na.rm = TRUE)
  rt_kurtosis <- kurtosis(behavior_data$RT, na.rm = TRUE)
  
  print(paste("Skewness:", round(rt_skewness, 4)))
  print(paste("Kurtosis:", round(rt_kurtosis, 4)))
  
  # Interpretation
  skew_interpretation <- if(abs(rt_skewness) < 1) {
    "Skewness is within acceptable range (|value| < 1)"
  } else if(rt_skewness > 1) {
    "Positive skew: distribution has a long right tail"
  } else {
    "Negative skew: distribution has a long left tail"
  }
  
  kurt_interpretation <- if(abs(rt_kurtosis - 3) < 2) {
    "Kurtosis is within acceptable range (|excess kurtosis| < 2)"
  } else if(rt_kurtosis > 5) {
    "Leptokurtic: distribution is more peaked with heavier tails than normal"
  } else {
    "Platykurtic: distribution is flatter with lighter tails than normal"
  }
  
  print(skew_interpretation)
  print(kurt_interpretation)
  
  # Normality tests
  print("\nNormality Tests:")
  
  # Shapiro-Wilk test (for n < 5000)
  if (length(behavior_data$RT) < 5000) {
    shapiro_test <- shapiro.test(na.omit(behavior_data$RT))
    print("Shapiro-Wilk test:")
    print(shapiro_test)
    print(paste("p-value:", shapiro_test$p.value))
    
    if (shapiro_test$p.value < 0.05) {
      print("Shapiro-Wilk test suggests non-normal distribution (p < 0.05)")
    } else {
      print("Shapiro-Wilk test suggests normal distribution (p >= 0.05)")
    }
  } else {
    print("Sample size too large for Shapiro-Wilk test, using Anderson-Darling instead")
  }
  
  # Anderson-Darling test (for larger samples)
  ad_test <- ad.test(na.omit(behavior_data$RT))
  print("\nAnderson-Darling test:")
  print(ad_test)
  print(paste("p-value:", ad_test$p.value))
  
  if (ad_test$p.value < 0.05) {
    print("Anderson-Darling test suggests non-normal distribution (p < 0.05)")
  } else {
    print("Anderson-Darling test suggests normal distribution (p >= 0.05)")
  }
  
  # Visual inspection
  print("\nCreating normality plots...")
  
  # Create QQ plot
  qq_plot_file <- file.path(result_path, "RT_QQ_plot.png")
  png(qq_plot_file, width = 800, height = 600)
  qqnorm(behavior_data$RT, main = "Q-Q Plot of RT Data")
  qqline(behavior_data$RT, col = "red")
  dev.off()
  print(paste("Q-Q Plot saved to:", qq_plot_file))
  
  # Create histogram
  hist_plot_file <- file.path(result_path, "RT_histogram.png")
  png(hist_plot_file, width = 800, height = 600)
  hist(behavior_data$RT, breaks = 30, main = "Histogram of RT Data", 
       xlab = "Reaction Time", prob = TRUE)
  curve(dnorm(x, mean = mean(behavior_data$RT, na.rm = TRUE), 
              sd = sd(behavior_data$RT, na.rm = TRUE)), 
        add = TRUE, col = "red", lwd = 2)
  dev.off()
  print(paste("Histogram saved to:", hist_plot_file))
  
  # Create density plot
  density_plot_file <- file.path(result_path, "RT_density_plot.png")
  png(density_plot_file, width = 800, height = 600)
  plot(density(behavior_data$RT, na.rm = TRUE), main = "Density Plot of RT Data")
  curve(dnorm(x, mean = mean(behavior_data$RT, na.rm = TRUE), 
              sd = sd(behavior_data$RT, na.rm = TRUE)), 
        add = TRUE, col = "red", lwd = 2)
  legend("topright", legend = c("RT Data", "Normal Distribution"), 
         col = c("black", "red"), lwd = 2)
  dev.off()
  print(paste("Density plot saved to:", density_plot_file))
  
  # Save normality analysis results
  normality_results <- data.frame(
    Statistic = c("Mean", "Median", "Standard Deviation", "Skewness", "Kurtosis",
                  "Shapiro-Wilk p-value", "Anderson-Darling p-value"),
    Value = c(
      mean(behavior_data$RT, na.rm = TRUE),
      median(behavior_data$RT, na.rm = TRUE),
      sd(behavior_data$RT, na.rm = TRUE),
      rt_skewness,
      rt_kurtosis,
      ifelse(length(behavior_data$RT) < 5000, shapiro_test$p.value, NA),
      ad_test$p.value
    ),
    Interpretation = c(
      "Central tendency",
      "Central tendency (robust to outliers)",
      "Spread of data",
      skew_interpretation,
      kurt_interpretation,
      ifelse(length(behavior_data$RT) < 5000, 
             ifelse(shapiro_test$p.value < 0.05, "Non-normal distribution", "Normal distribution"),
             "Sample size too large for test"),
      ifelse(ad_test$p.value < 0.05, "Non-normal distribution", "Normal distribution")
    )
  )
  
  # Save results to Excel
  normality_file <- file.path(result_path, "RT_normality_analysis.xlsx")
  write_xlsx(normality_results, normality_file)
  print(paste("\nNormality analysis results saved to:", normality_file))
  
  # Create Chinese conclusion text file
  chinese_conclusion_file <- file.path(result_path, "RT正态性分析结论.txt")
  
  # Prepare Chinese text interpretations
  chinese_skew <- if(abs(rt_skewness) < 1) {
    "偏度在可接受范围内 (|值| < 1)，分布呈现良好的对称性"
  } else if(rt_skewness > 1) {
    "呈现正偏度，分布右侧有长尾，数据分布向右偏"
  } else {
    "呈现负偏度，分布左侧有长尾，数据分布向左偏"
  }
  
  chinese_kurt <- if(abs(rt_kurtosis - 3) < 2) {
    "峰度在可接受范围内 (|超额峰度| < 2)，分布尖峭度正常"
  } else if(rt_kurtosis > 5) {
    "峰度过高，分布过于尖峭（厚尾），极端值较多"
  } else {
    "峰度过低，分布较为平坦（薄尾），极端值较少"
  }
  
  chinese_shapiro <- if(length(behavior_data$RT) < 5000) {
    if(shapiro_test$p.value < 0.05) {
      paste0("Shapiro-Wilk检验表明数据分布不符合正态分布 (p = ", round(shapiro_test$p.value, 4), " < 0.05)")
    } else {
      paste0("Shapiro-Wilk检验表明数据分布符合正态分布 (p = ", round(shapiro_test$p.value, 4), " >= 0.05)")
    }
  } else {
    "样本量过大，不适用于Shapiro-Wilk检验"
  }
  
  chinese_ad <- if(ad_test$p.value < 0.05) {
    paste0("Anderson-Darling检验表明数据分布不符合正态分布 (p = ", round(ad_test$p.value, 4), " < 0.05)")
  } else {
    paste0("Anderson-Darling检验表明数据分布符合正态分布 (p = ", round(ad_test$p.value, 4), " >= 0.05)")
  }
  
  chinese_overall <- if((abs(rt_skewness) < 1 && abs(rt_kurtosis - 3) < 2) || 
                      (length(behavior_data$RT) < 5000 && shapiro_test$p.value >= 0.05) || 
                      ad_test$p.value >= 0.05) {
    "综合多种检验方法，反应时（RT）数据近似符合正态分布。"
  } else {
    "综合多种检验方法，反应时（RT）数据偏离正态分布。"
  }
  
  chinese_transform <- if((abs(rt_skewness) >= 1 || abs(rt_kurtosis - 3) >= 2) || 
                         (length(behavior_data$RT) < 5000 && shapiro_test$p.value < 0.05) || 
                         ad_test$p.value < 0.05) {
    if(rt_skewness > 1) {
      "数据呈现正偏态，建议考虑对数变换（log transformation）、平方根变换（square root transformation）或倒数变换（inverse transformation）。"
    } else if(rt_skewness < -1) {
      "数据呈现负偏态，建议考虑平方变换（square transformation）或立方变换（cube transformation）。"
    } else {
      "考虑使用适当的数据变换使其更接近正态分布。"
    }
  } else {
    "数据近似服从正态分布，无需进行数据变换。"
  }
  
  # Write Chinese conclusion to text file
  cat(paste0(
    "反应时（RT）数据正态性分析结论\n",
    "==========================\n\n",
    "1. 基本统计量\n",
    "   - 样本量：", length(na.omit(behavior_data$RT)), "个\n",
    "   - 平均值：", round(mean(behavior_data$RT, na.rm = TRUE), 2), "\n",
    "   - 中位数：", round(median(behavior_data$RT, na.rm = TRUE), 2), "\n",
    "   - 标准差：", round(sd(behavior_data$RT, na.rm = TRUE), 2), "\n",
    "   - 偏度：", round(rt_skewness, 4), " (判据: |值| < 1 为可接受)\n",
    "   - 峰度：", round(rt_kurtosis, 4), " (判据: |值-3| < 2 为可接受)\n\n",
    "2. 正态性判断\n",
    "   - ", chinese_skew, "\n     偏度值为 ", round(rt_skewness, 4), "，判据: |值| < 1 为可接受\n",
    "   - ", chinese_kurt, "\n     峰度值为 ", round(rt_kurtosis, 4), "，判据: |值-3| < 2 为可接受\n",
    "   - ", chinese_shapiro, "\n     判据: p值 ≥ 0.05 表示符合正态分布\n",
    "   - ", chinese_ad, "\n     判据: p值 ≥ 0.05 表示符合正态分布\n\n",
    "3. 综合结论\n",
    "   ", chinese_overall, "\n",
    "   偏度判据: |", round(rt_skewness, 4), "| ", ifelse(abs(rt_skewness) < 1, "< 1 (满足)", ">= 1 (不满足)"), "\n",
    "   峰度判据: |", round(rt_kurtosis, 4), " - 3| = |", round(rt_kurtosis - 3, 4), "| ", 
    ifelse(abs(rt_kurtosis - 3) < 2, "< 2 (满足)", ">= 2 (不满足)"), "\n",
    "   Shapiro-Wilk检验: p = ", ifelse(length(behavior_data$RT) < 5000, round(shapiro_test$p.value, 4), "N/A"), 
    ifelse(length(behavior_data$RT) < 5000, 
           ifelse(shapiro_test$p.value >= 0.05, " ≥ 0.05 (满足)", " < 0.05 (不满足)"),
           " (样本量过大，不适用)"), "\n",
    "   Anderson-Darling检验: p = ", round(ad_test$p.value, 4),
    ifelse(ad_test$p.value >= 0.05, " ≥ 0.05 (满足)", " < 0.05 (不满足)"), "\n\n",
    "4. 数据变换建议\n",
    "   ", chinese_transform, "\n\n",
    "注意：正态性检验是对混合效应模型（Mixed Effects Model）的前提假设进行验证的重要步骤。\n",
    "如果数据严重偏离正态分布，建议采取适当的数据变换或考虑非参数方法。\n\n",
    "分析日期：", format(Sys.Date(), "%Y年%m月%d日"), "\n"
  ), file = chinese_conclusion_file, fileEncoding = "UTF-8")
  
  print(paste("\n中文结论文件已保存至:", chinese_conclusion_file))
  
  # Overall conclusion
  print("\nOverall conclusion about RT normality:")
  if ((abs(rt_skewness) < 1 && abs(rt_kurtosis - 3) < 2) || 
      (length(behavior_data$RT) < 5000 && shapiro_test$p.value >= 0.05) || 
      ad_test$p.value >= 0.05) {
    print("The RT data appears to be approximately normally distributed.")
  } else {
    print("The RT data appears to deviate from normal distribution.")
    
    # Suggestions for transformations if non-normal
    print("\nSuggested transformations for non-normal data:")
    if (rt_skewness > 1) {
      print("For positive skew, consider log transformation, square root, or inverse transformation.")
    } else if (rt_skewness < -1) {
      print("For negative skew, consider square or cube transformation.")
    }
  }
  
  # ============== Log-transformed RT Analysis ==============
  print("\n============= Log-transformed RT Normality Analysis =============")
  
  # Create log-transformed RT
  # Check for non-positive values before log transformation
  min_rt <- min(behavior_data$RT, na.rm = TRUE)
  
  if (min_rt <= 0) {
    print("Warning: RT contains zero or negative values. Adding a constant before log transformation.")
    # Add a constant to make all values positive
    log_rt <- log(behavior_data$RT + abs(min_rt) + 1)
  } else {
    log_rt <- log(behavior_data$RT)
  }
  
  # Add log-transformed RT to the dataframe
  behavior_data$log_RT <- log_rt
  
  print("Log transformation applied to RT data.")
  
  # Basic descriptive statistics for log-transformed RT
  log_rt_summary <- summary(behavior_data$log_RT)
  print("Log-transformed RT Summary Statistics:")
  print(log_rt_summary)
  
  # Skewness and Kurtosis for log-transformed RT
  log_rt_skewness <- skewness(behavior_data$log_RT, na.rm = TRUE)
  log_rt_kurtosis <- kurtosis(behavior_data$log_RT, na.rm = TRUE)
  
  print(paste("Log-transformed RT Skewness:", round(log_rt_skewness, 4)))
  print(paste("Log-transformed RT Kurtosis:", round(log_rt_kurtosis, 4)))
  
  # Interpretation for log-transformed RT
  log_skew_interpretation <- if(abs(log_rt_skewness) < 1) {
    "Skewness is within acceptable range (|value| < 1)"
  } else if(log_rt_skewness > 1) {
    "Positive skew: distribution has a long right tail"
  } else {
    "Negative skew: distribution has a long left tail"
  }
  
  log_kurt_interpretation <- if(abs(log_rt_kurtosis - 3) < 2) {
    "Kurtosis is within acceptable range (|excess kurtosis| < 2)"
  } else if(log_rt_kurtosis > 5) {
    "Leptokurtic: distribution is more peaked with heavier tails than normal"
  } else {
    "Platykurtic: distribution is flatter with lighter tails than normal"
  }
  
  print(log_skew_interpretation)
  print(log_kurt_interpretation)
  
  # Normality tests for log-transformed RT
  print("\nNormality Tests for Log-transformed RT:")
  
  # Shapiro-Wilk test for log-transformed RT (for n < 5000)
  if (length(behavior_data$log_RT) < 5000) {
    log_shapiro_test <- shapiro.test(na.omit(behavior_data$log_RT))
    print("Shapiro-Wilk test for log-transformed RT:")
    print(log_shapiro_test)
    print(paste("p-value:", log_shapiro_test$p.value))
    
    if (log_shapiro_test$p.value < 0.05) {
      print("Shapiro-Wilk test suggests non-normal distribution (p < 0.05)")
    } else {
      print("Shapiro-Wilk test suggests normal distribution (p >= 0.05)")
    }
  } else {
    print("Sample size too large for Shapiro-Wilk test, using Anderson-Darling instead")
  }
  
  # Anderson-Darling test for log-transformed RT
  log_ad_test <- ad.test(na.omit(behavior_data$log_RT))
  print("\nAnderson-Darling test for log-transformed RT:")
  print(log_ad_test)
  print(paste("p-value:", log_ad_test$p.value))
  
  if (log_ad_test$p.value < 0.05) {
    print("Anderson-Darling test suggests non-normal distribution (p < 0.05)")
  } else {
    print("Anderson-Darling test suggests normal distribution (p >= 0.05)")
  }
  
  # Visual inspection for log-transformed RT
  print("\nCreating normality plots for log-transformed RT...")
  
  # Create QQ plot for log-transformed RT
  log_qq_plot_file <- file.path(result_path, "log_RT_QQ_plot.png")
  png(log_qq_plot_file, width = 800, height = 600)
  qqnorm(behavior_data$log_RT, main = "Q-Q Plot of Log-transformed RT Data")
  qqline(behavior_data$log_RT, col = "red")
  dev.off()
  print(paste("Log-transformed RT Q-Q Plot saved to:", log_qq_plot_file))
  
  # Create histogram for log-transformed RT
  log_hist_plot_file <- file.path(result_path, "log_RT_histogram.png")
  png(log_hist_plot_file, width = 800, height = 600)
  hist(behavior_data$log_RT, breaks = 30, main = "Histogram of Log-transformed RT Data", 
       xlab = "Log(Reaction Time)", prob = TRUE)
  curve(dnorm(x, mean = mean(behavior_data$log_RT, na.rm = TRUE), 
              sd = sd(behavior_data$log_RT, na.rm = TRUE)), 
        add = TRUE, col = "red", lwd = 2)
  dev.off()
  print(paste("Log-transformed RT Histogram saved to:", log_hist_plot_file))
  
  # Create density plot for log-transformed RT
  log_density_plot_file <- file.path(result_path, "log_RT_density_plot.png")
  png(log_density_plot_file, width = 800, height = 600)
  plot(density(behavior_data$log_RT, na.rm = TRUE), main = "Density Plot of Log-transformed RT Data")
  curve(dnorm(x, mean = mean(behavior_data$log_RT, na.rm = TRUE), 
              sd = sd(behavior_data$log_RT, na.rm = TRUE)), 
        add = TRUE, col = "red", lwd = 2)
  legend("topright", legend = c("Log-transformed RT Data", "Normal Distribution"), 
         col = c("black", "red"), lwd = 2)
  dev.off()
  print(paste("Log-transformed RT Density plot saved to:", log_density_plot_file))
  
  # Save normality analysis results for log-transformed RT
  log_normality_results <- data.frame(
    Statistic = c("Mean", "Median", "Standard Deviation", "Skewness", "Kurtosis",
                  "Shapiro-Wilk p-value", "Anderson-Darling p-value"),
    Value = c(
      mean(behavior_data$log_RT, na.rm = TRUE),
      median(behavior_data$log_RT, na.rm = TRUE),
      sd(behavior_data$log_RT, na.rm = TRUE),
      log_rt_skewness,
      log_rt_kurtosis,
      ifelse(length(behavior_data$log_RT) < 5000, log_shapiro_test$p.value, NA),
      log_ad_test$p.value
    ),
    Interpretation = c(
      "Central tendency",
      "Central tendency (robust to outliers)",
      "Spread of data",
      log_skew_interpretation,
      log_kurt_interpretation,
      ifelse(length(behavior_data$log_RT) < 5000, 
             ifelse(log_shapiro_test$p.value < 0.05, "Non-normal distribution", "Normal distribution"),
             "Sample size too large for test"),
      ifelse(log_ad_test$p.value < 0.05, "Non-normal distribution", "Normal distribution")
    )
  )
  
  # Save results to Excel for log-transformed RT
  log_normality_file <- file.path(result_path, "log_RT_normality_analysis.xlsx")
  write_xlsx(log_normality_results, log_normality_file)
  print(paste("\nLog-transformed RT normality analysis results saved to:", log_normality_file))
  
  # Create Chinese conclusion text file for log-transformed RT
  log_chinese_conclusion_file <- file.path(result_path, "log_RT正态性分析结论.txt")
  
  # Prepare Chinese text interpretations for log-transformed RT
  log_chinese_skew <- if(abs(log_rt_skewness) < 1) {
    "偏度在可接受范围内 (|值| < 1)，分布呈现良好的对称性"
  } else if(log_rt_skewness > 1) {
    "呈现正偏度，分布右侧有长尾，数据分布向右偏"
  } else {
    "呈现负偏度，分布左侧有长尾，数据分布向左偏"
  }
  
  log_chinese_kurt <- if(abs(log_rt_kurtosis - 3) < 2) {
    "峰度在可接受范围内 (|超额峰度| < 2)，分布尖峭度正常"
  } else if(log_rt_kurtosis > 5) {
    "峰度过高，分布过于尖峭（厚尾），极端值较多"
  } else {
    "峰度过低，分布较为平坦（薄尾），极端值较少"
  }
  
  log_chinese_shapiro <- if(length(behavior_data$log_RT) < 5000) {
    if(log_shapiro_test$p.value < 0.05) {
      paste0("Shapiro-Wilk检验表明数据分布不符合正态分布 (p = ", round(log_shapiro_test$p.value, 4), " < 0.05)")
    } else {
      paste0("Shapiro-Wilk检验表明数据分布符合正态分布 (p = ", round(log_shapiro_test$p.value, 4), " >= 0.05)")
    }
  } else {
    "样本量过大，不适用于Shapiro-Wilk检验"
  }
  
  log_chinese_ad <- if(log_ad_test$p.value < 0.05) {
    paste0("Anderson-Darling检验表明数据分布不符合正态分布 (p = ", round(log_ad_test$p.value, 4), " < 0.05)")
  } else {
    paste0("Anderson-Darling检验表明数据分布符合正态分布 (p = ", round(log_ad_test$p.value, 4), " >= 0.05)")
  }
  
  log_chinese_overall <- if((abs(log_rt_skewness) < 1 && abs(log_rt_kurtosis - 3) < 2) || 
                          (length(behavior_data$log_RT) < 5000 && log_shapiro_test$p.value >= 0.05) || 
                          log_ad_test$p.value >= 0.05) {
    "综合多种检验方法，对数变换后的反应时（log_RT）数据近似符合正态分布。"
  } else {
    "综合多种检验方法，对数变换后的反应时（log_RT）数据仍偏离正态分布。"
  }
  
  log_chinese_transform <- if((abs(log_rt_skewness) >= 1 || abs(log_rt_kurtosis - 3) >= 2) || 
                             (length(behavior_data$log_RT) < 5000 && log_shapiro_test$p.value < 0.05) || 
                             log_ad_test$p.value < 0.05) {
    if(log_rt_skewness > 1) {
      "对数变换后数据仍呈现正偏态，可尝试其他变换方法如Box-Cox变换。"
    } else if(log_rt_skewness < -1) {
      "对数变换后数据呈现负偏态，可尝试其他变换方法。"
    } else {
      "考虑使用其他变换方法使数据更接近正态分布。"
    }
  } else {
    "对数变换有效改善了数据分布，数据现在近似服从正态分布。"
  }
  
  # Compare with original RT
  log_chinese_comparison <- if (
    (abs(log_rt_skewness) < abs(rt_skewness) && abs(log_rt_kurtosis - 3) < abs(rt_kurtosis - 3)) ||
    (length(behavior_data$RT) < 5000 && log_shapiro_test$p.value > shapiro_test$p.value) ||
    log_ad_test$p.value > ad_test$p.value
  ) {
    "对比原始反应时数据，对数变换明显改善了数据的正态性。"
  } else {
    "对比原始反应时数据，对数变换未能明显改善数据的正态性。"
  }
  
  # Write Chinese conclusion to text file for log-transformed RT
  cat(paste0(
    "对数变换后反应时（log_RT）数据正态性分析结论\n",
    "==========================\n\n",
    "1. 基本统计量\n",
    "   - 样本量：", length(na.omit(behavior_data$log_RT)), "个\n",
    "   - 平均值：", round(mean(behavior_data$log_RT, na.rm = TRUE), 2), "\n",
    "   - 中位数：", round(median(behavior_data$log_RT, na.rm = TRUE), 2), "\n",
    "   - 标准差：", round(sd(behavior_data$log_RT, na.rm = TRUE), 2), "\n",
    "   - 偏度：", round(log_rt_skewness, 4), " (判据: |值| < 1 为可接受)\n",
    "   - 峰度：", round(log_rt_kurtosis, 4), " (判据: |值-3| < 2 为可接受)\n\n",
    "2. 正态性判断\n",
    "   - ", log_chinese_skew, "\n     偏度值为 ", round(log_rt_skewness, 4), "，判据: |值| < 1 为可接受\n",
    "   - ", log_chinese_kurt, "\n     峰度值为 ", round(log_rt_kurtosis, 4), "，判据: |值-3| < 2 为可接受\n",
    "   - ", log_chinese_shapiro, "\n     判据: p值 ≥ 0.05 表示符合正态分布\n",
    "   - ", log_chinese_ad, "\n     判据: p值 ≥ 0.05 表示符合正态分布\n\n",
    "3. 综合结论\n",
    "   ", log_chinese_overall, "\n",
    "   偏度判据: |", round(log_rt_skewness, 4), "| ", ifelse(abs(log_rt_skewness) < 1, "< 1 (满足)", ">= 1 (不满足)"), "\n",
    "   峰度判据: |", round(log_rt_kurtosis, 4), " - 3| = |", round(log_rt_kurtosis - 3, 4), "| ", 
    ifelse(abs(log_rt_kurtosis - 3) < 2, "< 2 (满足)", ">= 2 (不满足)"), "\n",
    "   Shapiro-Wilk检验: p = ", ifelse(length(behavior_data$log_RT) < 5000, round(log_shapiro_test$p.value, 4), "N/A"), 
    ifelse(length(behavior_data$log_RT) < 5000, 
           ifelse(log_shapiro_test$p.value >= 0.05, " ≥ 0.05 (满足)", " < 0.05 (不满足)"),
           " (样本量过大，不适用)"), "\n",
    "   Anderson-Darling检验: p = ", round(log_ad_test$p.value, 4),
    ifelse(log_ad_test$p.value >= 0.05, " ≥ 0.05 (满足)", " < 0.05 (不满足)"), "\n",
    "   ", log_chinese_comparison, "\n\n",
    "4. 进一步建议\n",
    "   ", log_chinese_transform, "\n\n",
    "注意：对数变换是处理正偏态数据的常用方法，可以压缩数据的右侧尾部，使分布更接近正态。\n",
    "变换后的数据在统计分析中可能更符合正态性假设，但解释结果时需考虑变换的影响。\n\n",
    "分析日期：", format(Sys.Date(), "%Y年%m月%d日"), "\n"
  ), file = log_chinese_conclusion_file, fileEncoding = "UTF-8")
  
  print(paste("\n对数变换后反应时数据的中文结论文件已保存至:", log_chinese_conclusion_file))
  
  # Overall conclusion for log-transformed RT
  print("\nOverall conclusion about Log-transformed RT normality:")
  if ((abs(log_rt_skewness) < 1 && abs(log_rt_kurtosis - 3) < 2) || 
      (length(behavior_data$log_RT) < 5000 && log_shapiro_test$p.value >= 0.05) || 
      log_ad_test$p.value >= 0.05) {
    print("The log-transformed RT data appears to be approximately normally distributed.")
  } else {
    print("The log-transformed RT data still deviates from normal distribution.")
  }
  
  # Comparison with original RT
  if ((abs(log_rt_skewness) < abs(rt_skewness) && abs(log_rt_kurtosis - 3) < abs(rt_kurtosis - 3)) ||
      (length(behavior_data$RT) < 5000 && log_shapiro_test$p.value > shapiro_test$p.value) ||
      log_ad_test$p.value > ad_test$p.value) {
    print("Log transformation has improved the normality of the RT data.")
  } else {
    print("Log transformation did not substantially improve the normality of the RT data.")
  }
  
  # ============== Right Reaction Log-transformed RT Analysis ==============
  # 检查可能的reaction列
  reaction_columns <- names(behavior_data)[grep("react|response|resp|correct|accuracy", names(behavior_data), ignore.case = TRUE)]
  
  if (length(reaction_columns) > 0) {
    print("\n============= Right Reaction Log-transformed RT Normality Analysis =============")
    print("检测到可能的反应列:")
    print(reaction_columns)
    
    # 首先尝试精确名称"reaction"
    if ("reaction" %in% names(behavior_data)) {
      reaction_col <- "reaction"
    } else {
      # 否则使用找到的第一个可能列
      reaction_col <- reaction_columns[1]
    }
    
    print(paste("使用列:", reaction_col))
    
    # 查看该列的唯一值，帮助调试
    unique_values <- unique(behavior_data[[reaction_col]])
    print("该列的唯一值:")
    print(unique_values)
    
    # 以更健壮的方式筛选"right"反应 - 尝试多种可能的值
    right_data <- behavior_data %>% 
      filter(
        tolower(trimws(!!sym(reaction_col))) %in% 
          c("right", "correct", "1", "true", "yes", "t", "y", "r")
      )
    
    # 如果筛选结果为空，尝试数值型筛选
    if (nrow(right_data) == 0 && is.numeric(behavior_data[[reaction_col]])) {
      print("尝试将数值1视为正确反应...")
      right_data <- behavior_data %>% filter(!!sym(reaction_col) == 1)
    }
    
    # 如果筛选结果仍为空，再尝试使用其他列
    if (nrow(right_data) == 0 && length(reaction_columns) > 1) {
      for (alt_col in reaction_columns[-1]) {
        print(paste("尝试替代列:", alt_col))
        
        right_data <- behavior_data %>% 
          filter(
            tolower(trimws(!!sym(alt_col))) %in% 
              c("right", "correct", "1", "true", "yes", "t", "y", "r")
          )
        
        if (nrow(right_data) > 0) {
          reaction_col <- alt_col
          print(paste("成功使用替代列:", reaction_col))
          break
        }
        
        # 尝试数值型筛选
        if (is.numeric(behavior_data[[alt_col]])) {
          right_data <- behavior_data %>% filter(!!sym(alt_col) == 1)
          if (nrow(right_data) > 0) {
            reaction_col <- alt_col
            print(paste("成功使用替代列(数值型):", reaction_col))
            break
          }
        }
      }
    }
    
    print(paste("筛选后数据行数:", nrow(right_data)))
    
    # 检查筛选后是否有数据
    if (nrow(right_data) > 0 && "RT" %in% names(right_data)) {
      print(paste("成功筛选到", nrow(right_data), "行", reaction_col, "为正确反应的数据"))
      
      # 为right reaction数据创建log变换的RT
      # 检查负值或零值
      right_min_rt <- min(right_data$RT, na.rm = TRUE)
      
      if (right_min_rt <= 0) {
        print("警告: 正确反应RT包含零或负值。在log变换前添加常数。")
        # 添加常数使所有值为正
        right_log_rt <- log(right_data$RT + abs(right_min_rt) + 1)
      } else {
        right_log_rt <- log(right_data$RT)
      }
      
      # 将log变换的RT添加到数据框
      right_data$log_RT <- right_log_rt
      
      print("已对正确反应RT数据应用log变换。")
      
      # 基本描述性统计量
      right_log_rt_summary <- summary(right_data$log_RT)
      print("正确反应Log变换RT的统计摘要:")
      print(right_log_rt_summary)
      
      # 偏度和峰度
      right_log_rt_skewness <- skewness(right_data$log_RT, na.rm = TRUE)
      right_log_rt_kurtosis <- kurtosis(right_data$log_RT, na.rm = TRUE)
      
      print(paste("正确反应Log变换RT的偏度:", round(right_log_rt_skewness, 4)))
      print(paste("正确反应Log变换RT的峰度:", round(right_log_rt_kurtosis, 4)))
      
      # 解释
      right_log_skew_interpretation <- if(abs(right_log_rt_skewness) < 1) {
        "Skewness is within acceptable range (|value| < 1)"
      } else if(right_log_rt_skewness > 1) {
        "Positive skew: distribution has a long right tail"
      } else {
        "Negative skew: distribution has a long left tail"
      }
      
      right_log_kurt_interpretation <- if(abs(right_log_rt_kurtosis - 3) < 2) {
        "Kurtosis is within acceptable range (|excess kurtosis| < 2)"
      } else if(right_log_rt_kurtosis > 5) {
        "Leptokurtic: distribution is more peaked with heavier tails than normal"
      } else {
        "Platykurtic: distribution is flatter with lighter tails than normal"
      }
      
      print(right_log_skew_interpretation)
      print(right_log_kurt_interpretation)
      
      # 正态性检验
      print("\n正确反应Log变换RT的正态性检验:")
      
      # Shapiro-Wilk检验 (对于n < 5000)
      if (length(right_data$log_RT) < 5000) {
        right_log_shapiro_test <- shapiro.test(na.omit(right_data$log_RT))
        print("正确反应log变换RT的Shapiro-Wilk检验:")
        print(right_log_shapiro_test)
        print(paste("p值:", right_log_shapiro_test$p.value))
        
        if (right_log_shapiro_test$p.value < 0.05) {
          print("Shapiro-Wilk检验表明非正态分布 (p < 0.05)")
        } else {
          print("Shapiro-Wilk检验表明正态分布 (p >= 0.05)")
        }
      } else {
        print("样本量过大无法进行Shapiro-Wilk检验，改用Anderson-Darling检验")
      }
      
      # Anderson-Darling检验
      right_log_ad_test <- ad.test(na.omit(right_data$log_RT))
      print("\n正确反应log变换RT的Anderson-Darling检验:")
      print(right_log_ad_test)
      print(paste("p值:", right_log_ad_test$p.value))
      
      if (right_log_ad_test$p.value < 0.05) {
        print("Anderson-Darling检验表明非正态分布 (p < 0.05)")
      } else {
        print("Anderson-Darling检验表明正态分布 (p >= 0.05)")
      }
      
      # 可视化检验
      print("\n创建正确反应log变换RT的正态性图表...")
      
      # 创建QQ图
      right_log_qq_plot_file <- file.path(result_path, "right_log_RT_QQ_plot.png")
      png(right_log_qq_plot_file, width = 800, height = 600)
      qqnorm(right_data$log_RT, main = "正确反应Log变换RT数据的Q-Q图")
      qqline(right_data$log_RT, col = "red")
      dev.off()
      print(paste("正确反应log变换RT的Q-Q图已保存至:", right_log_qq_plot_file))
      
      # 创建直方图
      right_log_hist_plot_file <- file.path(result_path, "right_log_RT_histogram.png")
      png(right_log_hist_plot_file, width = 800, height = 600)
      hist(right_data$log_RT, breaks = 30, main = "正确反应Log变换RT数据的直方图", 
           xlab = "Log(反应时间)", prob = TRUE)
      curve(dnorm(x, mean = mean(right_data$log_RT, na.rm = TRUE), 
                  sd = sd(right_data$log_RT, na.rm = TRUE)), 
            add = TRUE, col = "red", lwd = 2)
      dev.off()
      print(paste("正确反应log变换RT的直方图已保存至:", right_log_hist_plot_file))
      
      # 创建密度图
      right_log_density_plot_file <- file.path(result_path, "right_log_RT_density_plot.png")
      png(right_log_density_plot_file, width = 800, height = 600)
      plot(density(right_data$log_RT, na.rm = TRUE), main = "正确反应Log变换RT数据的密度图")
      curve(dnorm(x, mean = mean(right_data$log_RT, na.rm = TRUE), 
                  sd = sd(right_data$log_RT, na.rm = TRUE)), 
            add = TRUE, col = "red", lwd = 2)
      legend("topright", legend = c("正确反应Log变换RT数据", "正态分布"), 
             col = c("black", "red"), lwd = 2)
      dev.off()
      print(paste("正确反应log变换RT的密度图已保存至:", right_log_density_plot_file))
      
      # 保存正态性分析结果
      right_log_normality_results <- data.frame(
        Statistic = c("Mean", "Median", "Standard Deviation", "Skewness", "Kurtosis",
                      "Shapiro-Wilk p-value", "Anderson-Darling p-value"),
        Value = c(
          mean(right_data$log_RT, na.rm = TRUE),
          median(right_data$log_RT, na.rm = TRUE),
          sd(right_data$log_RT, na.rm = TRUE),
          right_log_rt_skewness,
          right_log_rt_kurtosis,
          ifelse(length(right_data$log_RT) < 5000, right_log_shapiro_test$p.value, NA),
          right_log_ad_test$p.value
        ),
        Interpretation = c(
          "中心趋势",
          "中心趋势 (对离群值稳健)",
          "数据离散程度",
          right_log_skew_interpretation,
          right_log_kurt_interpretation,
          ifelse(length(right_data$log_RT) < 5000, 
                 ifelse(right_log_shapiro_test$p.value < 0.05, "非正态分布", "正态分布"),
                 "样本量过大无法检验"),
          ifelse(right_log_ad_test$p.value < 0.05, "非正态分布", "正态分布")
        )
      )
      
      # 保存结果到Excel
      right_log_normality_file <- file.path(result_path, "right_log_RT_normality_analysis.xlsx")
      write_xlsx(right_log_normality_results, right_log_normality_file)
      print(paste("\n正确反应log变换RT的正态性分析结果已保存至:", right_log_normality_file))
      
      # 创建中文结论文件
      right_log_chinese_conclusion_file <- file.path(result_path, "right_log_RT正态性分析结论.txt")
      
      # 准备中文文本解释
      right_log_chinese_skew <- if(abs(right_log_rt_skewness) < 1) {
        "偏度在可接受范围内 (|值| < 1)，分布呈现良好的对称性"
      } else if(right_log_rt_skewness > 1) {
        "呈现正偏度，分布右侧有长尾，数据分布向右偏"
      } else {
        "呈现负偏度，分布左侧有长尾，数据分布向左偏"
      }
      
      right_log_chinese_kurt <- if(abs(right_log_rt_kurtosis - 3) < 2) {
        "峰度在可接受范围内 (|超额峰度| < 2)，分布尖峭度正常"
      } else if(right_log_rt_kurtosis > 5) {
        "峰度过高，分布过于尖峭（厚尾），极端值较多"
      } else {
        "峰度过低，分布较为平坦（薄尾），极端值较少"
      }
      
      right_log_chinese_shapiro <- if(length(right_data$log_RT) < 5000) {
        if(right_log_shapiro_test$p.value < 0.05) {
          paste0("Shapiro-Wilk检验表明数据分布不符合正态分布 (p = ", round(right_log_shapiro_test$p.value, 4), " < 0.05)")
        } else {
          paste0("Shapiro-Wilk检验表明数据分布符合正态分布 (p = ", round(right_log_shapiro_test$p.value, 4), " >= 0.05)")
        }
      } else {
        "样本量过大，不适用于Shapiro-Wilk检验"
      }
      
      right_log_chinese_ad <- if(right_log_ad_test$p.value < 0.05) {
        paste0("Anderson-Darling检验表明数据分布不符合正态分布 (p = ", round(right_log_ad_test$p.value, 4), " < 0.05)")
      } else {
        paste0("Anderson-Darling检验表明数据分布符合正态分布 (p = ", round(right_log_ad_test$p.value, 4), " >= 0.05)")
      }
      
      right_log_chinese_overall <- if((abs(right_log_rt_skewness) < 1 && abs(right_log_rt_kurtosis - 3) < 2) || 
                                (length(right_data$log_RT) < 5000 && right_log_shapiro_test$p.value >= 0.05) || 
                                right_log_ad_test$p.value >= 0.05) {
        "综合多种检验方法，仅正确反应的对数变换反应时（log_RT）数据近似符合正态分布。"
      } else {
        "综合多种检验方法，仅正确反应的对数变换反应时（log_RT）数据仍偏离正态分布。"
      }
      
      right_log_chinese_transform <- if((abs(right_log_rt_skewness) >= 1 || abs(right_log_rt_kurtosis - 3) >= 2) || 
                               (length(right_data$log_RT) < 5000 && right_log_shapiro_test$p.value < 0.05) || 
                               right_log_ad_test$p.value < 0.05) {
        if(right_log_rt_skewness > 1) {
          "对数变换后数据仍呈现正偏态，可尝试其他变换方法如Box-Cox变换。"
        } else if(right_log_rt_skewness < -1) {
          "对数变换后数据呈现负偏态，可尝试其他变换方法。"
        } else {
          "考虑使用其他变换方法使数据更接近正态分布。"
        }
      } else {
        "对数变换有效改善了正确反应数据的分布，数据现在近似服从正态分布。"
      }
      
      # 写入中文结论到文本文件
      cat(paste0(
        "正确反应的对数变换反应时（log_RT）数据正态性分析结论\n",
        "==========================\n\n",
        "使用的反应列: ", reaction_col, "\n\n",
        "1. 基本统计量\n",
        "   - 样本量：", length(na.omit(right_data$log_RT)), "个\n",
        "   - 平均值：", round(mean(right_data$log_RT, na.rm = TRUE), 2), "\n",
        "   - 中位数：", round(median(right_data$log_RT, na.rm = TRUE), 2), "\n",
        "   - 标准差：", round(sd(right_data$log_RT, na.rm = TRUE), 2), "\n",
        "   - 偏度：", round(right_log_rt_skewness, 4), " (判据: |值| < 1 为可接受)\n",
        "   - 峰度：", round(right_log_rt_kurtosis, 4), " (判据: |值-3| < 2 为可接受)\n\n",
        "2. 正态性判断\n",
        "   - ", right_log_chinese_skew, "\n     偏度值为 ", round(right_log_rt_skewness, 4), "，判据: |值| < 1 为可接受\n",
        "   - ", right_log_chinese_kurt, "\n     峰度值为 ", round(right_log_rt_kurtosis, 4), "，判据: |值-3| < 2 为可接受\n",
        "   - ", right_log_chinese_shapiro, "\n     判据: p值 ≥ 0.05 表示符合正态分布\n",
        "   - ", right_log_chinese_ad, "\n     判据: p值 ≥ 0.05 表示符合正态分布\n\n",
        "3. 综合结论\n",
        "   ", right_log_chinese_overall, "\n",
        "   偏度判据: |", round(right_log_rt_skewness, 4), "| ", ifelse(abs(right_log_rt_skewness) < 1, "< 1 (满足)", ">= 1 (不满足)"), "\n",
        "   峰度判据: |", round(right_log_rt_kurtosis, 4), " - 3| = |", round(right_log_rt_kurtosis - 3, 4), "| ", 
        ifelse(abs(right_log_rt_kurtosis - 3) < 2, "< 2 (满足)", ">= 2 (不满足)"), "\n",
        "   Shapiro-Wilk检验: p = ", ifelse(length(right_data$log_RT) < 5000, round(right_log_shapiro_test$p.value, 4), "N/A"), 
        ifelse(length(right_data$log_RT) < 5000, 
               ifelse(right_log_shapiro_test$p.value >= 0.05, " ≥ 0.05 (满足)", " < 0.05 (不满足)"),
               " (样本量过大，不适用)"), "\n",
        "   Anderson-Darling检验: p = ", round(right_log_ad_test$p.value, 4),
        ifelse(right_log_ad_test$p.value >= 0.05, " ≥ 0.05 (满足)", " < 0.05 (不满足)"), "\n\n",
        "4. 进一步建议\n",
        "   ", right_log_chinese_transform, "\n\n",
        "注意：此分析仅针对正确反应的数据，共", nrow(right_data), "个观测值，占总数据的", 
        round(nrow(right_data)/nrow(behavior_data)*100, 2), "%。\n",
        "对数变换是处理正偏态数据的常用方法，可以压缩数据的右侧尾部，使分布更接近正态。\n",
        "变换后的数据在统计分析中可能更符合正态性假设，但解释结果时需考虑变换的影响。\n\n",
        "分析日期：", format(Sys.Date(), "%Y年%m月%d日"), "\n"
      ), file = right_log_chinese_conclusion_file, fileEncoding = "UTF-8")
      
      print(paste("\n正确反应对数变换反应时数据的中文结论文件已保存至:", right_log_chinese_conclusion_file))
      
      # 总体结论
      print("\n正确反应Log变换RT正态性的总体结论:")
      if ((abs(right_log_rt_skewness) < 1 && abs(right_log_rt_kurtosis - 3) < 2) || 
          (length(right_data$log_RT) < 5000 && right_log_shapiro_test$p.value >= 0.05) || 
          right_log_ad_test$p.value >= 0.05) {
        print("正确反应的log变换RT数据近似呈正态分布。")
      } else {
        print("正确反应的log变换RT数据仍偏离正态分布。")
      }
    } else {
      print("未找到正确反应的RT数据或筛选后无数据。")
      
      # 创建错误信息文件
      right_log_chinese_conclusion_file <- file.path(result_path, "right_log_RT正态性分析结论.txt")
      cat(paste0(
        "正确反应的对数变换反应时（log_RT）数据正态性分析结论\n",
        "==========================\n\n",
        "未能成功筛选到正确反应的数据。\n\n",
        "尝试使用的列: ", paste(reaction_columns, collapse=", "), "\n",
        "可能的原因:\n",
        "1. 数据中没有表示反应正确性的列\n",
        "2. 反应正确性的编码方式与预期不同\n",
        "3. 数据中没有正确反应的记录\n\n",
        "建议:\n",
        "1. 检查数据结构，确认正确反应的标记方式\n",
        "2. 查看各列的唯一值，确定哪一列包含反应正确性信息\n",
        "3. 修改代码中的筛选条件以匹配数据中的实际编码\n\n",
        "分析日期：", format(Sys.Date(), "%Y年%m月%d日"), "\n"
      ), file = right_log_chinese_conclusion_file, fileEncoding = "UTF-8")
      
      print(paste("\n正确反应数据筛选失败的错误信息已保存至:", right_log_chinese_conclusion_file))
    }
  } else {
    print("\n未在数据集中找到可能的反应列。")
    
    # 创建错误信息文件
    right_log_chinese_conclusion_file <- file.path(result_path, "right_log_RT正态性分析结论.txt")
    cat(paste0(
      "正确反应的对数变换反应时（log_RT）数据正态性分析结论\n",
      "==========================\n\n",
      "未在数据集中找到表示反应正确性的列。\n\n",
      "可能的原因:\n",
      "1. 数据结构中没有包含反应正确性的列\n",
      "2. 反应正确性列的命名方式与预期的模式不匹配\n\n",
      "建议:\n",
      "1. 检查数据结构，确认哪一列包含反应正确性信息\n",
      "2. 调整代码中搜索反应列的正则表达式以匹配实际列名\n\n",
      "数据集中的列名:\n",
      paste(names(behavior_data), collapse="\n"),
      "\n\n分析日期：", format(Sys.Date(), "%Y年%m月%d日"), "\n"
    ), file = right_log_chinese_conclusion_file, fileEncoding = "UTF-8")
    
    print(paste("\n未找到反应列的错误信息已保存至:", right_log_chinese_conclusion_file))
  }
} else {
  print("RT data not found in the dataset. Please check column names.")
  # Try to find possible RT columns with different naming
  possible_rt_cols <- names(behavior_data)[grep("rt|RT|reaction|time", names(behavior_data), ignore.case = TRUE)]
  if (length(possible_rt_cols) > 0) {
    print("Possible RT columns found:")
    print(possible_rt_cols)
  }
  
  # Create a message file if RT column not found
  chinese_conclusion_file <- file.path(result_path, "RT正态性分析结论.txt")
  cat(paste0(
    "反应时（RT）数据正态性分析结论\n",
    "==========================\n\n",
    "未在数据集中找到名为'RT'的列。请检查数据结构或列名。\n\n",
    "可能的反应时相关列：\n",
    paste(possible_rt_cols, collapse = "\n"),
    "\n\n分析日期：", format(Sys.Date(), "%Y年%m月%d日"), "\n"
  ), file = chinese_conclusion_file, fileEncoding = "UTF-8")
  
  print(paste("\n中文结论文件已保存至:", chinese_conclusion_file))
} 